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摘 要 : 本 文 基 于 二 阶 导数 的 四 阶 Pad6 型 紧 致 差分 逼近 式 ， 并 结合 原 方程 本 身 ， 得 到 了 三 
维 Helmholtz 方 程 的 一 种 四 阶 精 度 的 隐 式 紧 致 差分 格式 ， 该 格式 在 每 个 空间 方向 上 只 涉及 到 
三 个 点 处 的 未 知 量 及 其 二 阶 导 数值 。 边 界 处 对 于 二 阶 导数 的 离散 格式 利用 四 阶 显 式 偏心 格式 。 然 
后 ， 利 用 Richardson 外 推 法 、 算 子 插值 法 及 二 阶 导 数 在 边界 点 处 的 六 阶 显 式 偏心 格式 ， 将 本 文 
构造 的 格式 精度 提高 到 六 阶 。 最 后 ， 通 过 数值 实验 验证 了 本 文 方法 的 精确 性 和 可 靠 性 。 
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1 引言 


Helmholtz 方 程 在 许多 工程 领域 都 会 遇 到 ， 如 电磁 场 中 的 波导 问题 、 薄 膜 振 动 问题 以 及 
海洋 工程 中 水 波 衍 射 问 题 等 都 是 由 Helmholtz 方 程控 制 的 。Helmholtz 方 程 的 数值 解法 有 多 
种 ， 常 见 的 有 有 限 单元 法 和 有 限 差 分 法 及 其 各 种 迭代 方法 Q3， 其 中 有 限 单元 法 应 用 最 为 
普遍 ， 但 是 ， 有 限 元 法 的 计算 精度 随 着 方程 波 数 的 增 大 而 急剧 降低 。 而 传统 的 差分 方法 精 
度 较 低 ， 要 想得到 精确 的 结果 ， 则 必须 将 网 格 精细 章 分 ， 如 此 是 很 不 经 济 的 ， 于 是 研究 者 
倾向 于 发 展 各 种 高 精度 的 差分 方法 8-6l 。Manohar 和 Stephenson?! 提出 了 二 维 Hemholtz 方程 
的 一 种 高 精度 差分 格式 ，Singer Al Turkell4] 提出 了 关于 三 维 Hemholtz 方 程 的 一 种 高 精度 紧 
致 差分 格式 ; Steijl 和 Hoeijmakers”! 基于 Laplace 算 子 的 高 阶 紧 致 离散 形式 ， 推 导 得 到 了 三 
维 Helmholtz 方 程 的 一 种 高 阶 紧 致 差分 格式 。 王 前 东 等 中 给 出 了 变 系数 Helmholtz 方 程 的 六 阶 
精度 校正 法 ， 此 校正 过 程 几乎 不 增加 工作 量 ， 而 校正 后 精度 比 未 校正 提高 了 四 阶 。 

本 文采 用 二 阶 导数 的 四 阶 Pad6 型 紧 致 差分 逼近 式 ， 结 合 原 方程 本 身 , 得 到 了 三 维 Helmholtz 
方程 的 一 种 四 阶 精 度 的 隐 式 紧 致 差分 格式 ， 该 格式 在 每 个 空间 方向 上 只 涉及 到 三 个 点 处 的 未 知 
量 及 其 二 阶 导数 值 ， 边 界 处 对 于 二 阶 导数 的 离散 采用 文献 中 提出 的 四 阶 显 式 偏心 格式 。 然 后 
基于 Richardson 外 推 法 、 算 子 插值 法 及 二 阶 导数 在 边界 点 处 的 六 阶 显 式 偏心 格式 ， 将 本 文 构 造 
的 格式 精度 整体 提高 到 六 阶 。 
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2 高 阶 隐 式 紧 致 差分 格式 
考虑 如 下 的 三 维 Helmholtz 方 程 


| Au(z,y,2) + 2u(z,y, 2) = f(z,9,2), (2,9,2) €D, 





(1) 
u(x, y, 2) = g(x,y, z), (x,y,z) € OD, 
其 中 人 是 三 维 Laplace 算 子 ， 刀 是 求解 区 域 ，09D 是 区 域 边界 ，#k 是 为 常数 。 现 取 计 算 区 域 
为 [0, lz] x [0, ly] x [0, 1z]， 每 个 空间 方向 取 等 步 长 。 记 
le o oL 

he = N7’ hy = he = N 

Nz, Ny, Nz 为 网 格 等 分 数 ， 且 均 为 偶数 ，hs, hy, he 分别 为 z,y,z 方 向 的 网 格 步 长 。 
根据 二 阶 导 数 的 四 阶 Padé 型 紧 致 差分 格式 四 


1 5 1 Ui4t1 jk — 2ZUi ik + Uili 
T(r) + g (Uaa) + 75 Uk = a + O(hs), (2) 
zx 
1 5 1 Uij+l,k — 2Ui,jk F Ui j-1,k 
Ta Man istik + g (Uyy)iak + Toy) = = a E Olh), (3) 
y 
5 1 Ui j,k+1 — Qui, ik 十 Ui j,k- 
Tg (U2) ktl + g (Uez)iak + Tz (Uez)ijk-1 = E + O(hs), (4) 
zZz 
以 及 
Ui j,k = Ui,j,k: (5) 


将 (2) x 8 十 (3) x £4(4) x 94(5) x e2， 并 利用 原 方程 (1) 式 ， 略 去 高 阶 项 后 整理 可 得 
pg (uit +u- Lat) + gpg (tistre + Mes Lk) 


12 12 12 


Ui j,k+1 + Ui, j,k- 1)+ (r? T A r Sa ) vik 


ey 
+n 5h2 5h2 5h2 


1 
= a5 [(Ure)i+1,j,k + (Ua)i bik + (Uyy)ij+i,k + (Uyy)i,j—1,k 


+ (tez)ig,b+1 + (Uzz)ijk-1] + Seg e- (6) 


方程 (6) 就 是 求解 三 维 Helmholtz 方程 的 高 阶 隐 式 紧 致 差分 格式 ， 记 为 HOIC。 由 推导 过 程 
可 知 其 截断 误差 为 DO( 奏 十 网 十 彤 )， 其 中 关于 未 知 量 的 二 阶 导数 可 由 (2), (3) 和 (4) 式 给 出 ， 由 
于 其 值 在 边界 点 处 是 未 知 的 ， 因 此 ， 为 了 保持 格式 的 整体 精度 达到 四 阶 ， 利 用 文献 如 给 出 其 
人 如 二 阶 导 数 在 z 方 向 左右 边界 点 处 的 离散 格式 为 


(tzz)0,j,4 = 225uo jik 一 770u1 j,k + 1070w2,;, 天 一 780u3 j,k + 305u4 jik 一 50us dy k) (7) 


sore ( 


(Ura) N。 (225un,,,j,k = 770UN, —1,j,k + 1070uN, —2,7,k 780UN, —3,5,k 


jk = o 
+ 305UN,—4,j,k — 50uN,—5,j,k) (8) 
二 阶 导数 在 y 方 向 和 z 方 向 上 边界 点 处 的 离散 格式 可 类 似 给 出 。 
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3 Richardson 外 推 法 


本 文 的 计算 空间 包括 三 部 分 ， BAD, yh.’ Dam oho Diaby he? Pha hon BO 
FAA HE ARTI Ze 9 Ohi + hy + ha) BE, Din, on, on, ATL RTI EW O((2he)* + 
(2hy)* + (2hz)*) BRZI, DP, PAA UTR Ze Wy Oh, + hs + hs) 的 解 空间 。 


对 四 阶 紧 致 差分 格式 ，Richardson 外 推 法 的 计算 公式 为 


16u}, 5.9, — uch 
GP yp = EE OS 15, RS N/2. (9) 


为 了 书写 方便 ， 用 DERE Dhan’ Din RE Dh on, on.» DIRE DS nao BRE 
说 在 实际 计算 中 并 不 要 求 网 格 是 等 距 的 ， 可 以 是 非 等 距 网 格 。 具 体 算法 描述 如 下 : 


1) MAA ub p E Di, uh, © Di Mib, EDS, JAMRA eps; 


2) 利用 初 值 几 jx € Di, 求 细 网 格 上 未 知 量 及 其 二 阶 导 数 的 四 阶 精度 解 , BD (usr), (uyy), 
(uzz)” 和 ub € Di. 
(a) RÆ (2), (7) 和 (8) 三 式 中 系数 和 矩阵 为 (Ns — 1) x (Ne — 1) 的 三 对 角 线 性 方程 组 ， 求 
得 (uze)” E€ Di; 
(b) 同样 方法 求 得 (wyy)* € D4 和 (wss)* € D4; 
(c) 通过 方程 (6) 可 得 (10) 式 ， 用 于 求解 D4 上 z 方 向 wr 的 值 ， 只 需求 解 (10) 式 xz 方向 上 
形成 的 (Ns - 1) x (Ns 一 1) 的 三 对 角 线 型 方程 组 


12 12 12 
2 h 
sya (ut Page HUR rit) + (4 ~ Bad ~ Bre ~ ana) ia 


= h h h 
= 5 [Curr) tisk + (Uzz)i Lk + (Uyy)ij+i,k + (uyy) jik + (urz) ee + (zz) kar] 


6 h h 
+ fij,k 一 5h2 (ww Ui j+l,k + ut 9-1 = sia (ul 十 ths Ped) 


1=1,2,.…., (Nz — 1), j=1,2,.…,(N, ~—1), k= 1,2,...,(N; — 1). (10) 


3) PADRE u, © D 包 可 以 求 得 粗 网 格 上 未 知 量 及 其 二 阶 导 数 的 四 阶 精度 解 ， 
即 (t4e2)?”, (uyy), (Uz2) Alu?’ © D4,. WR 2); 


4) 利用 (urs)”, (uyy)”, (uze)? € Dt 和 (uss), (uyy)? ”, (uz2)™ € Då, Kükk e€ D6 求 
(ae) (ip) (ŭz2)” = D$. 

(a) 利用 (Uae )3i,25,2k € Dh, (ura)? i € D3, 和 Richardson 外 推 , 求 得 (tha )3i,25,2k € Dh: 

(b) 利用 (tira )3;.93,0% € DÊ, 他 JE Di 和 (2) 式 ， 由 算 子 插值 法 ， 得 (thw )3;—1,25,2% EDR: 


(c) 利用 ttt 5 € Di, (ten )b,0j—-1,0k € D$, (tas) N, 2j—1,2k E€ D8$ 和 (2) 式 ， 由 算 子 插值 
法 ， fF (üre) 9; 19k € D$, 其 中 (izz)b 2;_1,2k € D$ 和 (tine) N, 2;_1,2k ED}; 采用 文献 [7] 对 
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二 阶 导 数 的 六 阶 逼 近 式 。 


a 1 7469 223 879 949 
(iizz)o,jk = m (= gg Wk Tq Wak + Gy jk — Tg “3k 


201 1019 7 
+ 4lua jk 一 Jo 57k + Taq 67 一 TAE il: (11) 
: 469 223 879 949 
(tex) Na, j,k = m ae UNa jik 一 g UNe—lik + gg UNe—2ik — Tg UNe 3k 
201 1019 7 
+ 41UN,-4,5,k — Fp UNz-57 + eq “Nek ~ TQ UNe—Tik k): (12) 


(a) 利用 
Ui j,k E Df, (iaz)072k-1 E Dj, (Han) N, 426-1 € Dh, 

MOR MURE (de) jana € DË, SE Gea)h jana € DEM ee), yarı € DE RAK 
BK [7] St SRS Wr SK (11) 和 (12) 式 ; 

(e) 同样 方法 求 得 (tow) jk E D$ Al (人 ez) kE DE: 

5) Aub, E€ Dh, utp € Dan, (iizz)”, (tyy)” Aah, € DI RARE a” 的 六 阶 精度 

解 ， 即 g? Ee D$, 

(a) 利用 uhi 2j,2k € Di M uzh e € Pow 应 用 Richardson SME» RI üh, 2;2k € Dh 

(b) 利用 (tex), (diyy)”, (za) iP yk E D$, ag, ,2j,2k © Dj 和 (13) 式 ， 应 用 算 子 插值 
法 ， 求 得 雄 1,27,2k € Dj: 


Eog z : 
Baa +j + ayn) + (« ~ BRR BRB BAB” 
= ii u ul i a i 
= 10 (CaS eee + (ize) ijk + (好 十 (Gy) oak + (za) Pj kt + (zz) tj er) 


6 6 
+fijk — = +u U; + Wij, 
tJ Baa “i t,jt1,k tj- 1k) 一 aa (ut t,j,k+1 es 1) 
k = 2,4,--- ,(Nz — 2), j =2,4,--- ,(Ny — 2), 4=1,3,-:-,(Nz—1); (13) 


(c) 利用 (tes) (ŭyy)”, (ŭz2)", Thik € Df, Thi 2j,2k € Dĵ, Ühi—1,2j,2k € D$ 和 (6) 式 ， 
同样 方法 求 得 识 。;_1,2 E€ DBs 
(d) 利用 (iizs)*， (fyy)”, (tize)*, ühik € Df, ab, ,2j,2k © D$, ùh- ,2j,2k © D$, ùl Koy ak E 
Ds 和 (6) 式 ， 同 样 方法 求 得 ;2k_1 © Dh 
6) ”收敛 性 判断 : 如 果 ic > eps， 即 入 两 次 迭代 的 误差 大 于 等 于 eps， 转 到 步 2); 
7) 输出 结果 ， 停 止 程序 运 


4 数值 实验 


在 立方 体 空间 0 < zx,y,z < 1 内 ， 考 察 三 维 Helmboltz 方 程 如 下 两 个 有 精确 解 的 问题 。 本 文 
将 构造 的 基于 Richardson 外 推 法 的 Helmholtz 方 程 的 高 阶 隐 式 紧 致 差分 格式 记 为 RHOIC。 

问题 1 u(z,y,z) = sin(z)sin(2y)sin(3z)。 

问题 2 ulz, y, z) =g yz entot, 
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表 1 至 表 4 分 别 给 出 了 k= 1.6, 上 = 3.2 和 网 格 等 分 数 N = 16, 32, 64 时 ， 问 题 1 和 问题 2 采 
用 中 心 差分 格式 (CDS)、 四 阶 紧 致 格式 (FOC)S! 及 本 文 格式 RHOIC 数值 计算 结果 的 最 大 误 
差 |ell。、 均 方 根 误差 ||ell> 以 及 收敛 阶 rate 的 比较 ， 其 中 FOC 格式 是 文献 3] 二 维 方法 向 三 维 
的 直接 推广 。 

由 计算 结果 可 以 看 出 ， 对 于 两 个 问题 ， 本 文 格式 RHOIC 达到 了 六 阶 精度 ， 对 问题 1 和 问 
题 2， 本 文 格式 RHOIC 的 精确 性 高 于 FOC 格 式 ， 并 且 本 文 格式 RHOIC 的 精度 明显 高 于 中 心 
差分 格式 。 昌 然 本 文 没 有 列 出 CDS 格 式 h = 1/256 时 的 数值 解 ， 但 是 本 文 格式 RHOIC 在 粗 
网 格 h = 1/16 下 的 精确 度 仍 高 于 CDS 格 式 在 精细 网 格 h = 1/256 下 的 结果 。 例 如 对 问题 1， 
当 = 1.6 时 ， 本 文 格式 RHOIC 在 粗 网 格 h = 1/16 下 的 最 大 误差 为 6.81 x 10-6, mi CDS Hat 
在 精细 网 格 h = 1/256 下 的 最 大 误差 为 1.38 x 10-5. 


表 1: 问题 1 三 种 格式 在 k = 1.6 时 的 误差 和 收敛 阶 的 比较 














FOC 格式 加 
lello llell2 
1.56(-5)  5.62(-6) 
9.90(-7) 3.50(-7) 3.98 
6.81(-8) 2.91(-8) 4.00 


CDS 格 式 
llelleo llelle 
3.53(-3) ”1.25(-3) 
8.80(-4) 3.11(-4) 2.00 
2.20(-4) 7.77(-5) 2.00 


RHOIC 格式 
































llelloo llell2 
6.81(-6)  4.45(-6) 

7.63(-8)  8.11(-8) 6.48 
7.71(-10) 1.32(-9) 6.63 








表 2: 问题 1 三 种 格式 在 k& = 3.2 时 的 误差 和 收敛 阶 的 比较 













FOCHA! 
llelloo llelle rate 
2.22(-5) 7.84(-6) 
1.38(-6) 4.89(-7) 4.00 
8.63(-8)  3.05(-8) 4.00 


CDS 格 式 
llelloo llell2 

16 | 4.93(-3) 1.74(-3) 
32 | 1.23(-3)  4.35(-3) 2.00 
64 | 3.07(-4)  1.09(-4) 2.00 







RHOIC 格式 
llelloo llelle rate 
7.34(-6) 6.07(-6) 
7.78(-8) 1.03(-8) 6.56 
7.82(-10)  1.61(-10) 6.64 























表 3: 问题 2 三 种 格式 在 k = 1.6 时 的 误差 和 收敛 阶 的 比较 



















FOC 格式 四 
llelloo llell2 
7.47(-6)  2.43(-6) 
4.68(-7) 1.52(-7) 4.00 
2.93(-8) 9.53(-9) 4.00 






CDS 格 式 
llellw llell2 
1.18(-3)  3.27(-4) 
32 | 3.05(-4)  8.29(-4) 1.95 
64 | 7.65(-5) 2.08(-5) 1.99 


RHOIC 格式 


llelle rate 












llell~ 
6.39(-8)  6.84(-8) 

1.53(-9) 3.21(-9) 5.38 
3.21(-11) 1.12(-10) 5.61 
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表 4: 问题 2 三 种 格式 在 = 3.2 时 的 误差 和 收 敏 阶 的 比较 


llelloo llell2 rate | llellw llell2 
1.40(-3)  4.22(-4) 9.21(-6)  3.22(-6) 
3.54(-4) 1.07(-4) 1.98 | 5.77(-7) 2.02(-7) 4.00 
8.87(-5)  2.68(-5) 3.61(-8) 1.26(-8) 
















RHOIC 格式 




























6.88(-8) 9.90(-8) 
1.52(-9) 4.22(-9) 5.50 
2.92(-11) 1.51(-10) 
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High-order Implicit Compact Difference Scheme for Solving the 
Three-dimensional Helmholtz Equation 


GE Yong-bin', LIU Guo-tao? 


(1- Institute of Applied Mathematics and Mechanics, Ningxia University, Yinchuan 750021; 
2- The 46st Institute of the sixth Academy of CASIC, Hohhot 010010) 

Abstract: Based on the Padé scheme of second-order partial derivatives and combined with the 
original differential equation, a fourth-order implicit compact difference scheme is proposed for solving 
the three-dimensional Helmholtz equation. Only three points and their second-order derivative values 
are needed on each spacial direction. The fourth-order explicit difference schemes are used to construct 
the same order discretization of boundary conditions. Then, the accuracy of the fourth-order implicit 
compact difference scheme is upgraded to sixth-order by using the Richardson extrapolation technique 
and operator interpolation scheme. The sixth-order explicit difference schemes of second-order partial 
derivatives on the boundaries are also used. Finally, numerical experiments are given to prove the 
efficiency and reliability of the present method. 

Keywords: Helmholta equation; high accuracy; implicit; compact difference scheme; Richardson 
extrapolation 
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